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Abstract. On the basis of additive schemes (splitting schemes) we construct 
efficient numerical algorithms to solve approximately the initial-boundary value 
problems for systems of time-dependent partial differential equations (PDEs). 
In many applied problems the individual components of the vector of unknowns 
are coupled together and then splitting schemes are applied in order to get a 
simple problem for evaluating components at a new time level. Typically, 
the additive operator-difference schemes for systems of evolutionary equations 
are constructed for operators coupled in space. In this paper we investigate 
more general problems where coupling of derivatives in time for components 
of the solution vector takes place. Splitting schemes are developed using an 
additive representation for both the primary operator of the problem and the 
operator at the time derivative. Splitting schemes are based on a triangular 
two-component representation of the operators. 



Introduction 

Mathematical modeling of applied problems is based on approximate solving 
boundary value problems for systems of time-dependent PDEs. To construct nu- 
merical algorithms for solving such problems, we develop approximations to the 
equations taking into account the corresponding initial and boundary conditions. 
The approximation in space is conducted using finite difference schemes, finite el- 
ement procedures or finite volume methods [3 [321 EI]. Special attention should 
be given to approximations in time for numerical solving problems with systems 
of equations [HEJH]. ^ n Edition to the general conditions of approximation and 
stability we should keep in mind the questions of computational implementation 
of the constructed schemes. To solve the corresponding grid problem at the new 
time level in an efficient way, special additive operator-difference schemes (splitting 
schemes) are in common use [51 IT51II?!] . 

In approximate solving initial-boundary value problems for multidimensional 
PDEs, a transition to a chain of simpler problems allows us to construct economi- 
cal difference schemes — we have splitting with respect to spatial variables. In some 
cases it is useful to split into the subproblems of different nature — we have split- 
ting into physical processes [SHU]. Recently, active discussions are concerned with 
regionally additive schemes (domain decomposition schemes), which are oriented 
to parallel computers jlO[fT5]. 

The additive schemes for vector problems can be treated as a separate class. 
The schemes of this type can be used to develop efficient numerical algorithms for 
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solving time-dependent systems of PDEs. A typical situation is the case where the 
individual components of the unknown vector are interconnected and it is difficult 
to derive a simple problem for evaluating vector components at the new time level. 

Various classes of additive schemes are developed for vector problems (see, e.g. 
[23]). For parabolic and hyperbolic systems of equations with a self-adjoint ellip- 
tic operator the locally one-dimensional additive schemes are constructed in |13j 
using the principle of regularization for difference schemes. To construct efficient 
splitting schemes for the systems of equations, the Samarskii alternating triangular 
method can be employed, which is generally treated as an iterative method [T3"I[Trj] . 
This approach is implemented, in particular, in |21j for the dynamic problems of 
elasticity, and in [22j — for the problems of an incompressible fluid with a variable 
viscosity. Additive schemes for time-dependent vector equations of first and second 
order are constructed in [19j for problems of electrodynamics. 

Traditionally, additive operator-difference schemes for systems of evolutionary 
equations are in common use to construct approximations to equations with opera- 
tors interconnected in space. But in some cases there does exist a coupling between 
the time derivatives of the individual components of the solution vector. There- 
fore, it is necessary to design additive operator-difference schemes with splitting an 
operator at the time derivative. Theory and practice of constructing the splitting 
schemes of such kind is right at the beginning of the development. In fact, the first 
work on the splitting schemes for the problems with an additive representation of 
an operator at the time derivative for the evolutionary equations of first order is 
the paper [20] ■ New vector additive schemes are proposed and investigated in it 
using splitting the operator at the time derivative into a sum of positive definite 
self-adjoint operators. Unfortunately, these schemes cannot be directly applied to 
the systems of evolutionary equations with a coupling between the time derivatives. 

In this paper we construct splitting schemes using an additive representation for 
both the primary operator of the problem and the operator at the time derivative. 
The schemes are based on a triangular two-component representation of the opera- 
tors and are applied to systems of evolutionary equations. The work is organized as 
follows. In Section 1 we formulate the Cauchy problem for a system of evolutionary 
equations of first order. We also provide an example of the equation system with 
PDEs, which can be related to the considered Cauchy problem with the correspond- 
ing approximations in space. The standard two-level operator-difference scheme is 
discussed in Section 2. Section 3 deals with the construction of additive schemes 
via the triangular splitting of the primary operator of the problem. The general 
problem with splitting both the primary operator of the problem and the operator 
at the time derivative is presented in Section 4. 



Let H a , a = 1,2, ...,p be finite-dimensional real Hilbert (Euclidean) spaces 
where the scalar product and the norm are (•, -) a and || • \\ a , a = 1, 2, . . . ,p, re- 
spectively. The individual components of the solution are denoted by u a (t), a — 
1,2, ... ,p for any t (0 < t < T, T > 0). We search the solution for the system of 
evolutionary equations of first order: 



1. Problem formulation 



(1.1) 




a = 1,2,..., p. 



/3=1 
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Here f a {t) £ £2(0, T; H a ), a — 1,2, ...,f> are specified, and B a p, A a p are the 
linear constant (independent of t) operators acting from Hp into H a (A a p : Hp — > 
H a , B a p : Hp — > H a ) for all a — 1,2,..., p. The system of equations (jl.ip is 
supplemented with the initial conditions 

(1.2) Ua (0)=v° a , a = 1,2,..., p. 

We treat the system of equations (jl.lj) as a single evolutionary equation for the 
vector u = {m, u-i, . . . , u p }: 

(1.3) B^+Au = f(i), 0<<<T, 

dt 

where f = {fi, f 21 ■■■■> fp}> ancl f° r the elements of the operator matrices A and B 
we have the representation 

A = {A a p}, B = {B a p}, a, = 1,2,..., p. 

On the direct sum of spaces [6] H = Hi © H 2 © • • • © H p we set 

p p 

( U > V ) = ^2(Ua,V a ) a , ||uf = 

Taking into account f|1.2jl . we define 

(1.4) u(0) = v°, 

where v° = {v®, v®, ■ ■ ■ , v®}. 

We consider the Cauchy problem (|1.3[) . (11.41) . under the condition that the op- 
erators A and B are self-adjoint and positive definite in H: 

(1.5) A = A* > SaE, S a >0, B = B* > 5 B B, 5 B > 0, 

where E is the identity operator in H. The self-adjointness is associated with the 
fulfillment of 

A a p = A*p a , B a p = Bp a , a,/3 = 1,2, ...,p 

for operators of the original system of equations (|1.1[) . 

Here is the simplest a priori estimate for the solution of the Cauchy problem 
(|1.3[) . (|1.4p . We will use it as a guide in investigating the corresponding operator- 
difference schemes. For D = D* > we use the notation Hd for a space H with 
the scalar product (y,w)d = (Dy,w) and norm ||y ||n> = (Dyjy) 1 / 2 . 

Scalar ly multiplying both sides of equation (11.31) in H by — , we obtain 

dt 

„du du\ Id., . / du 
B — , — + Au, u) = f, — 

dt' dt ) 2dV ' V dt 

Taking into account (|1.5[) and 

f)^(Bf.§^^-*.fl 

we derive the inequality 

d„ ,., 1 



B-i- 



We get from it the following a priori estimate 

(1-6) l|u(t)[[i<||v°||i + ijr t ||f(6l)|||_ lt », 
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which expresses the stability of the solution of problem (|1.3|) , (|1.4p with respect to 
the initial data and right-hand side. 

Systems of evolutionary equations, similar to result from the approxima- 

tion in space to many applied problems. Let us consider here some typical examples 
do not specifying the corresponding spaces of continuous and grid functions. 

First, we mention the system of coupled parabolic equations of second order 
describing mass transfer in multicomponent media [4}ll8j. The solution is sought 
in a bounded domain Q, u Q (x, t), x G f2: 

P du P P 

^^q/}(x)-^ - ^div(/c Q/3 (x)gradw /3 ) + ^ r a p(*)up = / a (x,f), 

0=1 0=1 p=i 

xefl, 0<t<T, a = 1,2,..., p. 

The coefficients r a p are associated with the reaction processes whereas k a g describe 
the diffusion phenomena (main-term diffusion at a = /3 and cross-term diffusion 
coefficients at a ^ 0). For the multicomponent media we have b a/ 3 = S a pb a , where 
5 a p is the Kronecker symbol. 

The second example concerns a fluid motion in the porous media. The governing 
equations for a flow in fractured porous media employ the multiple porosity model 
(see, for example, [2,3 ). In this case u a {x.,t) is the dynamic pore pressure in an 
p-porosity model. For these models it is principal that b a p ^ 5 a pb a and r a p =/= 
(the Barcnblatt model). 

2. Scheme with weights 

To solve approximately operator-differential problem (|1.3|) . (|1.4p . we use the 
standard scheme with weights |13j . Introduce a uniform grid in time 

57 T = cj t U {T} = {t n = nr, n = 0, 1, ...,N, tN = T} 

and suppose y" = y(t n ), t n = tit. Let us approximate equation (|1.3I) via the 
following two-level difference scheme 

(2.1) B y " +1 ~ y " + A(ay" +1 + (1 - a)y") = y", 

r 

where a is a numerical parameter (weight) within < a < 1 and, for example, 
(p n = f(crt" +1 + (1 — u)t n ). For simplicity, we restrict ourselves to the case of the 
equal weight for all equations of systems . In view of (|1.4p we supplement (|2.1[) 
with the initial condition 

(2.2) y° = v°. 

A detailed study of the schemes with weights (the necessary and sufficient conditions 
for stability, the choice of a norm) was conducted in [TJirrS]. Here we restrict 
ourselves to the simplest estimate of stability for operator-difference scheme (|2.1[) . 
(|2.2[) . Estimate (|1.6p will serve us as a guide. 

Theorem 2.1. If a > 1/2, then operator- difference scheme \2.1\) is absolutely 
stable in Ha and for the difference solution the following level-wise estimate is 
valid: 

(2.3) iiy" +1 iii^iiy"Hi + p^iiJ B+(CT _, )rArl . 
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Proof. We write scheme (|2.1j) in the form 



/ l\ \ v™ +1 - v™ v n+1 +v™ 
B+ ( t r--)rA)3L_JL + A 3L_t3L s=v n. 

Scalar ly multiplying both sides of this equation in H by 2(y™ +1 — y"), we obtain 
the equality 

2.[[B + L-lW y " +1 - y " y " +1 



2) } T T 

/ v n+l 

(Ay" +1 ,y" +1 ) - (Ay",y") = 2r <p«, Z 



Using the inequality 

n+l 




we get the required estimate (|2.3[) . □ 



Estimate (|2.3[) is just a grid analog of estimate (|1.6I) and provides the uncondi- 
tional stability of the difference scheme with weights (|2.1[) , (|2.2p under the natural 
conditions cr > 1/2. Considering the corresponding problem for the error, we obtain 
convergence of the solution of operator-difference problem (|2.1[) , (|2.2[) to the solu- 
tion of operator-differential problem (TOf . (fOJ) in H A at cr > 1/2 C((2ct-1)t+t 2 ). 
If a = 1/2, we have the second-order convergence with respect to r. 

Operator-difference scheme (|2.ip can be written in the canonical form of the 
two-level schemes [13] 

n+l _ n 

(2.4) (B + cttA)^ £_ + Ay" = <p". 

r 

The transition to a new time level requires to solve the problem 

(B + crrA)y n+1 = rf> n . 

For the original problem (|1.1[) , (jl.2|) we have to solve the system of coupled equations 

p 

J2( B <*0 + °TA aP ) y ; +1 =C, a = 1,2, . . .,p. 

13=1 



Various iterative methods can be used for it [1T ( [T6 ] . 

Another opportunity here is to take into account the specific features of the 
considered time-dependent problems and to construct splitting schemes, where the 
transition to a new time level is based on the solution of more simpler problems. 
For problems of type (jl.ip . (|1.2|) it is natural to employ the splitting schemes, where 
the transition to the new time level is performed via solving the problems 

(B aa +o-TA aa )yl +1 a= 1,2, ...,p. 

This means that we have to invert only the diagonal part of the operator matrix 
B + crrA in our computations. 
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3. Schemes with a diagonal operator B 



We start from the case where the problem of inversion of the operator B does not 
exist. Such a situation occurs when the operator matrix B at the time derivatives 
is diagonal: 



(3.1) 



B a p = 5 a pB ai a = 1,2, • • • ,p. 



This class of problems appears in simulation of mass transfer in multicomponent 
media. In this case the components of the solution vector are coupled due to the 
elements A a », a ^ (3 of the operator matrix A. 

Following to [r3J [2TJ [52] , let us construct the additive operator-difference 
schemes using the triangular splitting of the operator A: 



(3.2) 



Ai 



In the additive representation (|3.2[) we have 
Ai = 



(\A X1 





• \ 


A 2 i 


±^22 •• 












\ A pi 


A p2 ■■ 





A 2 = 





A 12 ■ 


• M P \ 





\A 22 ■ 


■ A 2p 






Ap—Xp 


V o 





2^PP } 



Instead of (|2.4p we use the scheme 

,,n+l 



(3.3) 



B^ 



- Ay" = V « 



where the operator B has the following factorized form 
(3.4) B = (B + o-rAi)B- 1 (B + <TrA 2 ). 



Scheme (|3.3p . (J3T4J) is an operator- matrix analog of the Samarskii alternating tri- 
angle method [T3]- 

Taking into account (|3.2p . due to self-adjointness and positive defmiteness of B 
at a > 1/2 we have 



B = B + cttA + o- 2 r 2 AiB~ 1 A 2 , B = B* > B 
Similar to Theorem 12. 11 we prove the following statement. 



<ttA. 



Theorem 3.1. If a > 1/2, then factorized operator- difference scheme Il3.ty) - pn$\ ) 
is absolutely stable in Ha and for the difference solution the following level-wise 
estimate is valid: 



T n+1 |2 



< 



- 2 w 



(B + (cr-i)TA+cr 2 T 2 A 1 B- 1 A 2 ) 1 ' 



The computational implementation of scheme (13.2I) - (|3.4[) for solving problem 
p.ip . (|1.2p . (13.11) can be conducted using the sequence of simpler problems: 



B a 



a-A aa )y^ 2 = r o - 



o--A aa )yl +1 =r a , a = 1,2,. 



. ,p. 



As the scheme with weights (|2.ip . factorized scheme (|3.2|l - (|3.4|l has the convergence 
of second order at a — 1/2 and of first order at any other weights. 
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4. Generalizations 

If problem (|1.3j) . (|1.4j) has a non-diagonal operator B (B a p ^ 5 a /3B a ), then ad- 
ditive operator-difference schemes can be constructed using the triangular splitting 
for both the operator A and the operator B. Similarly (|3.2I) . we assume that 

(4.1) B = B 1 +B 2 , B*=B 2 . 

Taking into account (|3.2j) . (|4.ip . we write the scheme with weights (|2.4I) as follows 

n+l _ n 

(4.2) C± y — + Ay" = <p n , 

T 

where 

(4.3) C = C 1 +C 2 , Ci=Bi+<7tAi, C 2 = B 2 + cttA 2 . 
In view of (( 1 . 5 j) we have 

C 1= C 2 , C a > ^(S B + <jtS a )E. 

The operator Ci + C 2 can be represented like this 

C a + C 2 = i(d + eE)(C 2 + eE) - i(C x - eE)(C 2 - eE) 
2e 2e 

at any e > 0. The value of e will be defined a little bit later. Instead of two-level 

scheme (|4.2j) we employ the three-level scheme 

1 v " +1 - v" 

(4.4) (c 1+e E)(C 2 +EE)i y - 

Is t 

1 v " - v"^ 1 
-(C a - eE)(C 2 - eE)* y - + Ay" = <p n . 

2e t 

The primary potential advantage of this scheme in compare with scheme (|4.2j) is 
that its implementation is based on the inversion of the factorized operator (Ci + 
eE)(C 2 + eE) at the new time level. 
Taking into account that 

y«+i _ y« y"^ 1 — y" — 1 t y"^ 1 — 2y" + y"~ 1 

t = Yt + 2 V 2 ' 

y" _ y™ _1 y"^ 1 — y" _1 t y™^ 1 — 2y" + y" 1 



t 2r 
we can write scheme (14.41) as follows 



v n+l _ v n— 1 v"^ 1 — 2v" 4- v" -1 

(4.5) C± — i + D± V— i + Ay" = <p'\ 

2r r z 

where 

D = ^(CiC 2 +e 2 E). 

2e 

In view of C = B + err A we verify directly that operator-difference scheme 
approximates equation (|1 .3[) with the first-order accuracy with respect to r at all 
e = 0(1). We now formulate the sufficient conditions of stability for this scheme. 
A comprehensive study on the stability of three-level schemes with self-adjoint 
operators was done in [13TH5] . We do not use here the general results on the 
stability of operator-difference schemes from the above works. Similar to [20], we 
obtain immediately the simplest estimates of stability with respect to the initial 
data and right-hand side. 
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Taking into account 



y" = \(y n+1 + 2y" + y"" 1 ) - \(y n+1 - 2y n + y"- 1 ) 



we write (14.51) as 



v n+l _ v n— 1 / _2 \ v n+l _ o n _l Ar ri — 1 

(4.6) C- - + [D- — A^ y 



2r V 4 

n+l _ 2v n i v n-l 

4 

Introducing 

v" = ^(y" + y"- 1 ), w" = ^ i — 

we can rewrite (14.61) in the form 



w «+i _ w « i 

(4.7) C- + R + -A(v" +1 + y") = V n , 

2 t 2 

where 

r 2 

R = D A. 

4 

Scalarly multiplying both sides of (|4.7[) by 

2(v" +1 - v") = r(w n+1 + w"), 

we get the equality 

(4.8) ^(C(w" +1 + w n ),w n+1 +w") + (R(w™ +1 - w"),w™ +1 +w n )- 

(A(v™ +1 + v"),v" +1 - v n ) = V,w" +1 +w"). 
For the right-hand side we use the estimate 

( t p",w™+ 1 +w n ) < i(C(w n+1 +w") + i(C-V,¥' n )- 
This yields from (|4.8[) the inequality 

(4-9) £ n+1 <5„ + I(C-V n ,«P n ), 

where we use the notation 

£„ = (Av",v") + (Rw",w"). 



Inequality (|4.9p will be the desired a priori estimate, if we show that £ n defines the 
squared norm of the difference solution. Due to the positivity of A, it is sufficient 
to require the non-negativity of the operator R. 
With the above-mentioned notations we have 

2 

R = -p(Ci + eE)(C 2 + eE) + -^(Ci - eE)(C 2 - eE) - ^A > 

4e 4e 4 

2 2 

-^(C x + eE)(C 2 + eE) - ^A > T —{a - 1)A. 

4e 4 4 

Thus, R > at a > 1. The result of our considerations is the following statement. 
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Theorem 4.1. If a > 1, then operator- difference scheme L3.2\) . j4-l\ )> U-4\ ) * s 

absolutely stable, and for the difference solution a priori estimate holds with 



f 



R 



The proven estimate (|4.9[) ensures the stability of operator-difference scheme 
()3.2|) . (|4.1|) . (|4.4|) with respect to the initial data and right-hand side. It can be 
treated as a more complex analog of estimate (|2.3|) and it agrees with estimate (|1.6j) 
for the solution of the problem (jl.3p . (I1.4[) . 

The computational implementation of scheme (|4.4I) involves the solution of the 
grid problem 

(4.10) (C 1 +£E)(C 2 + £ E)y n+1 = t/>" 

at the new time level. Introducing the auxiliary unknown y n+1 / 2 , we have for (|4.10|) 
the following representation 

(C 1 +eE)y"+ 1 / 2 = -i/A 

(C 2 +eE)y" +1 =y ,l+1/2 . 

Taking into account the triangular structure of the operators (Ci and (C2, we 
sequentially solve the problem 

(B aa + arA aa + 2eE a )yl+ 1 ' 2 = fa 

{B aa + arA aa + 2eE a )yl +1 = fa a = l,2,...,p, 
where E a is the identity operator in H a . 
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